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Abstract 



We describe algorithms for finding the regression of t, a sequence of values, to the closest sequence 
s by mean squared error, so that s is always increasing (isotonicity) and so the values of two consecu- 
tive points do not increase by too much (Lipschitz). The isotonicity constraint can be replaced with a 
unimodular constraint, where there is exactly one local maximum in s. These algorithm are generalized 
from sequences of values to trees of values. For each scenario we describe near-linear time algorithms. 
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1 Introduction 



Let M be a triangulation of a polygonal region P C R 2 in which each vertex is associated with a real 
valued height (or elevation). Linear interpolation of vertex heights in the interior of each triangle of M 
defines a piecewise-linear function t : P —> R, called a height function. A height function (or its graph) 
is widely used to model a two-dimensional surface in numerous applications (e.g. modeling the terrain of 
a geographical area). With recent advances in sensing and mapping technologies, such models are being 
generated at an unprecedentedly large scale. These models are then used to analyze the surface and to 
compute various geometric and topological properties of the surface. For example, researchers in GIS are 
interested in extracting river networks or computing visibility or shortest-path maps on terrains modeled as 
height functions. These structures depend heavily on the topology of the level sets of the surface and in 
particular on the topological relationship between its critical points (maxima, minima, and saddle points). 
Because of various factors such as measurement or sampling errors or the nature of the sampled surface, 
there is often plenty of noise in these surface models which introduces spurious critical points. This in 
turn leads to misleading or undesirable artifacts in the computed structures, e.g., artificial breaks in river 
networks. These difficulties have motivated extensive work on topological simplification and noise removal 
through modification of the height function t into another one s : M — > R that has the desired set of critical 
points and that is as close to t as possible (7J|2TJ|26l|25j[33l. A popular approach is to decompose the surface 
into pieces and modify each piece so that it has a unique minimum or maximum [27]. In some applications, 
it is also desirable to impose the additional constraint that the function s is Lipschitz; see below for further 
discussion. 

Problem statement. Let M = (V, A) be a planar graph with vertex set V and arc (edge) set A C V x V. 
We may treat M as undirected in which case we take the pairs (u, v) and (v, u) as both representing the 
same undirected edge connecting u and v. Let 7 > be a real parameter. A function s : V — > R is called 

(L) ^-Lipschitz if (it, v) € A implies s(v) — s(u) < 7. 

Note that if M is undirected, then Lipschitz constraint on an edge (u, v) £ A implies \s(u) — s(v)\ < 7. For 
an undirected planar graph M = (V, A), a function s : V —> R is called 

(U) unimodal if s has a unique local maximum, i.e. only one vertex v € V such that s(v) > s(u) for all 
(u,v) £ A. 

For a directed planar graph M = (V, A), a function s : V — > R is called 

(I) isotonic if (it, v) € A implies s(it) < s(i>)Q 

For an arbitrary function t : V — > R and a parameter 7, the ^-Lipschitz unimodal regression (j-LUR) of t is a 
function s : V — > R that is 7-Lipschitz and unimodal on M and minimizes \\s — 1\\ 2 = ^2 v& y{s(v) —t(v)) 2 . 
Similarly, if M is a directed planar graph, then s is the ^-Lipschitz isotonic regression (j-LIR) of t if s 
satisfies (L) and (I) and minimizes \\s — t\\ 2 . The more commonly studied isotonic regression (IR) and 
unimodal regression E SI [23j are the special cases of LIR and LUR, respectively, in which 7 = 00, and 
therefore only the condition (I) or (U) is enforced. 

Given a planar graph M, a parameter 7, and t : V — > R, the LIR (resp. LUR) problem is to compute 
the 7-LIR (resp. 7-LUR) of t. In this paper we propose near-linear-time algorithms for the LIR and LUR 
problems for two special cases: when M is a path or a tree. We study the special case where M is a path 
prior to the more general case where it is tree because of the difference in running time and because doing 
so simplifies the exposition to the more general case. 

1 A function s satisfying the isotonicity constraint (I) must assign the same value to all the vertices of a directed cycle of M 
(and indeed to all vertices in the same strongly connected component). Therefore, without loss of generality, we assume M to be a 
directed acyclic graph (DAG). 
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Related work. As mentioned above, there is extensive work on simplifying the topology of a height 
function while preserving the geometry as much as possible. Two widely used approaches in GIS are the 
so-called flooding and carving techniques fl] QT1 |25j |26]|. The former technique raises the height of the 
vertices in "shallow pits" to simulate the effect of flooding, while the latter lowers the value of the height 
function along a path connecting two pits so that the values along the path vary monotonically. As a result, 
one pit drains to the other and thus one of the minima ceases to exist. Various methods based on Laplacian 
smoothing have been proposed in the geometric modeling community to remove unnecessary critical points; 
see [5J [7J QjO |2T| and references therein. For example, Ni et al. |2T1 proposed the so-called Morse-fairing 
technique, which solves a relaxed form of Laplace's equation, to remove undesired critical points. 

A prominent line of research on topological simplification was initiated by Edelsbrunner et al. lfl4l[T3l 
who introduced the notion of persistence; see also lfl2l l35l l34l . Roughly speaking, each homology class 
of the contours in sublevel sets of a height function is characterized by two critical points at one of whom 
the class is born and at the other it is destroyed. The persistence of this class is then the height difference 
between these two critical points and can be thought of as the life span of that class. The persistence 
of a class effectively suggests the "significance" of its defining critical points. Efficient algorithms have 
been developed for computing the persistence associated to the critical points of a height function and for 
simplifying topology based on persistence tfTBl 1711. Edelsbrunner et al. iflBI and Attali et al. proposed 
algorithms for optimally eliminating all critical points of persistence below a threshold where the error is 
measured as 1 1 s — £ 1 1 oo = max„ e y \ s(v ) — t(v)\. No efficient algorithm is known to minimize [|s — t\\2. 

The isotonic -regression (IR) problem has been studied in statistics H |6l [23j since the 1950s. It has 
many applications ranging from statistics l|29l to bioinformatics 0, and from operations research EUl to 
differential optimization iTTTl . Ayer et. al. H famously solves the IR problem on paths in 0{n) time using 
the pool adjacent violator algorithm (PAVA). The algorithm works by initially treating each vertex as a 
level set and merging consecutive level sets that are out of order. This algorithm is correct regardless of the 
order of the merges G4l . Brunk fl9] and Thompson [ 32 ] initiated the study of the IR problem on general 
DAGs and trees, respectively. Motivated by the problem of finding an optimal insurance rate structure over 
given risk classes for which a desired pattern of tariffs can be specified, Jewel |[T9l introduced the problem 
of Lipschitz isotonic regression on DAGs and showed connections between this problem and network flow 
with quadratic cost function^] Stout PTll solves the UR problem on paths in 0(n) time. Pardalos and 
Xue [22] give an 0(n log n) algorithm for the IR problem on trees. For the special case when the tree is 
a star they give an 0(n) algorithm. Spouge et al. |[28l give an 0(n 4 ) time algorithm for the IR problem 
on DAGs. The problems can be solved under the L\ and norms on paths [31] and DAGs as well, 
with an additional logn factor for L\. Recently Stout [30] has presented an efficient algorithm for istonotic 
regression in a planar DAG under the norm. To our knowledge no polynomial-time algorithm is known 
for the UR problem on planar graphs, and there is no prior attempt on achieving efficient algorithms for 
Lipschitz isotonic/unimodal regressions in the literature. 



Our results. Although the LUR problem for planar graphs remains elusive, we present efficient exact al- 
gorithms for LIR and LUR problems on two special cases of planar graphs: paths and trees. In particular, we 
present an 0(n log n) algorithm for computing the LIR on a path of length n (Section|4]), and an 0(n log n) 
algorithm on a tree with n nodes (Section 6|). We present an 0(n log 2 n) algorithm for computing the LUR 
problem on a path of length n (Section [5 1. Our algorithm can be extended to solve the LUR problem on 
an unrooted tree in 0(nlog 3 n) time (Section [7}. The LUR algorithm for a tree is particularly interesting 



2 It may at first seem that Jewel's formulation of the LIR s of an input function t on the vertices of a DAG is more general than 
ours in that he requires that for each e = (u, v) £ A, s(v) £ [s(u) — \ e , s(u) + j e ] where \ e ,y e £ R + are defined separately for 
each edge (in our formulation A e = and 7 e = 7 for every edge e). Moreover, instead of minimizing the L2 distance between s and 
t he requires w v (s(v) — t(v)) 2 to be minimized where w v S R + is a constant assigned to vertex v. However, all our algorithms 
in this paper can be modified in a straight-forward manner to handle this formulation without any change in running-times. 
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because of its application in the aforementioned carving technique ltTTll26ll27l . The carving technique mod- 
ifies the height function along a number of trees embedded on the terrain where the heights of the vertices 
of each tree are to be changed to vary monotonically towards a chosen "root" for that tree. In other words, 
to perform the carving, we need to solve the IR problem on each tree. The downside of doing so is that the 
optimal IR solution happens to be a step function along each path toward the root of the tree with potentially 
large jumps. Enforcing the Lipschitz condition prevents sharp jumps in function value and thus provides a 
more natural solution to the problem. 

Section|3]presents a data structure, called affine composition tree (ACT), for maintaining a xy-monotone 
polygonal chain, which can be regarded as the graph of a monotone piecewise-linear function F : R — > R. 
Besides being crucial for our algorithms, ACT is interesting in its own right. A special kind of binary 
search tree, an ACT supports a number of operations to query and update the chain, each taking O(logn) 
time. Besides the classical insertion, deletion, and query (computing F(x) or F~ 1 (x) for a given x G R), 
one can apply an Interval operation that modifies a contiguous subchain provided that the chain remains 
x-monotone after the transformation, i.e., it remains the graph of a monotone function. 



2 Energy Functions 

On a discrete set U, a real valued function s : U —> R can be viewed as a point in the \U\ -dimensional 
Euclidean space in which coordinates are indexed by the elements of U and the component s u of s associated 
to an element u G U is s(u). We use the notation R 17 to represent the set of all real- valued functions defined 
on U. 

Let M = (V, A) be a directed acyclic graph on which we wish to compute 7-Lipschitz isotonic regres- 
sion of an input function t G M. v . For any set of vertices U C V, let M.[U] denote the subgraph of M induced 
by U, i.e. the graph (U,A[U]), where A[U] = {(u,v) G A : u,v G U}. The set of 7-Lipschitz isotonic 
functions on the subgraph M[£7] of M constitutes a convex subset of M. u , denoted by T(M, U). It is the 
common intersection of all half-spaces determined by the isotonicity and Lipschitz constraints associated 
with the edges in A[U], i.e., s u < s v and s v < s u + 7 for all (u, v) G A[U]. 

For U C V we define E v : R v -> R as 

= £(*(«) -W- 

Thus, the 7-Lipschitz isotonic regression of the input function t is a = arg min sgr ( M y\ Ey(s). For a 
subset [/CV and v G U define the function : R — > R as 

E MU]( X ) = min E u(s). 

Lemma 2.1 For any U C V and v G U, the function E^^ is continuous and strictly convex. 
Proof: Given x, y G R and < a < 1, consider the functions 

z x = arg min Ejj(z) 
zer(M,U);z(v)=x 

z v = arg min E\j{z) 
zer(M,U);z(y)=y 

z a = az x + (1 - a)z y . 
The strong convexity of Ejj implies that 

aEu(z x ) + (1 - a)Eu(zy) > Eu{z a ). 
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Figure 1: Left: the graph of a monotone piecewise linear function F (in solid black). Vertices are marked by hollow 
circles. The dashed curve is the result of applying the linear transform ipo(q) = Mq where M = ( y 3 y^)- The 

gray curve, a translation of the dashed curve under vector c = ( z\ )> is the result of applying Affine(^>) to F where 
ip(q) = Mq + c. Right: Interval(?/;, t~ , t + ), only the vertices of curve F whose x-coordinates are in the marked 
interval [r~, r + ] are transformed. The resulting curve is the thick gray curve. 

Furthermore, the convexity of T(M, U) implies that z a is also in T(M, U). Since by definition z a (v) = 
ax + (1 — a)y, we deduce that 

Eu(z a ) > min Ejj(z). 

zGF(M,u);z(v)=ax+(l-a)y 

Thus, the minimum is precisely E^(ax + (1 — a)y). This proves that is strictly convex on E and 
therefore continuous (and in fact differentiable except in countably many points). ■ 

3 Affine Composition Tree 

In this section we introduce a data structure, called affine composition tree (ACT), for representing an xy- 
monotone polygonal chain in IR 2 , which is being deformed dynamically. Such a chain can be regarded as the 
graph of a piecewise-linear monotone function F : E — > E, and thus is bijective. A breakpoint of F is the 
^-coordinate of a vertex of the graph of F (a vertex of F for short), i.e., a b G E at which the left and right 
derivatives of F disagree. The number of breakpoints of F will be denoted by \F\. A continuous piecewise- 
linear function F with breakpoints b\ < • ■ ■ < b n can be characterized by its vertices qi = (pi, F(bi)),i = 
1, . . . , n together with the slopes /x_ and fi + of its left and right unbounded pieces, respectively, extending 
to — oo and +oo. An affine transformation of M 2 is a map <fi : M? — ► M 2 , q t— ► M • g + c where M is a 
nonsingular 2x2 matrix, (a linear transformation) and c £ I 2 is a translation vector — in our notation we 
treat q G M 2 as a column vector. 

An ACT supports the following operations on a monotone piecewise-linear function F with vertices 
qi = (bi,F(bi)),i = 1, ... ,n: 

1. Evaluate(o) and Evaluate -1 ^): Given any a G E, return F(o) or F _1 (a). 

2. iNSERT(g): Given a point g = (x,y) G M 2 insert g as a new vertex of F. If x G this 
operation removes the segment qiqi+i from the graph of F and replaces it with two segments qiq and 
qqi+i, thus making x a new breakpoint of F with F(x) = y. If x < b\ or x > then the affected 
unbounded piece of F is replaced with one parallel to it but ending at q and a segment connecting q 
and the appropriate vertex of F and /x_ remain intact). We assume that F(6j) < y < F(bi+i). 
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3. Delete(6): Given a breakpoint b of F, removes the vertex (b,F(b)) from F; a delete operation 
modifies F in a manner similar to insert. 

4. Affine(^): Given an affine transformation tp : M. 2 —>■ M 2 , modify the function F to one whose graph 
is the result of application of ip to the graph of F. See Figure [TjLeft). 

5. iNTERVALp^, t - , r + ): Given an affine transformation V : K. 2 — ^ ^ 2 an d r _ ,r + G R and p G 
{x,y}, this operation applies ip to all vertices v of F whose p-coordinate is in the range [t~,t + ]. 
Note that Affine(^) is equivalent to Interval p (V>, — oo, +oo) for p = {x, y}. See Figure [TjRight). 

It must be noted that Affine and Interval are applied with appropriate choice of transformation 
parameters so that the resulting chain remains xy-monotone. 

An ACT 7 = 7(F) is a red-black tree that stores the vertices of F in the sorted order, i.e., the zth node 
of 7 is associated with the zth vertex of F. However, instead of storing the actual coordinates of a vertex, 
each node z of T stores an affine transformation <p z : q M z -q + c z . If zq, . . . , z^ = zis the path in 7 from 
the root zq to z, then let <& z (q) = 4> Z0 ((f> Zl (• • . 4>z k (q) ■■■))■ Notice that <& z is also an affine transformation. 
The actual coordinates of the vertex associated with z are (q x , q y ) = 3> z (0) where = (0, 0). 

Given a value b G R andp G {x, y}, let Pred p (6) (resp. SiJCC p (6)) denote the rightmost (resp. leftmost) 
vertex q of F such that the p-coordinate of q is at most (resp. least) b. Using ACT T, Pred p (6) and SUCC P (6) 
can be computed in 0(log n) time by following a path in T, composing the affine transformations along the 
path, evaluating the result at 0, and comparing its p-coordinate with b. We can answer EVALUATE(a) 
queries by determining the vertices q~ = PRED x (a) and q + = Succ^a) of F immediately preceding 
and succeeding a and interpolating F linearly between q~ and q + ; if a < b± (resp. a > bk), then F(a) 
is calculated using F(b±) and (resp. F(bk) and fj, + ). Since 6 _ and b + can be computed in O(logn) 
time and the interpolation takes constant time, Evaluate(o) is answered in time O(logra). Similarly, 
Evaluate -1 (a) is answered using Pred j/ (o) and Succ y (a). 

A key observation of ACT is that a standard rotation on any edge of 7 can be performed in O(l) time by 
modifying the stored affine transformations in a constant number of nodes (see Figure [2]) based on the fact 
that an affine transformation : q i— > M -q + c has an inverse affine transformation 4>~ 1 : q i— > M~ 1 • (q — c) ; 
provided that the matrix M is invertible. A point q G M. 2 is inserted into T by first computing the affme 
transformation $ u for the node u that will be the parent of the leaf z storing q. To determine <p z we solve, 
in constant time, the system of (two) linear equations $u(0 z (O)) = <7- The result is the translation vector 
c z . The linear transformation M z can be chosen to be an arbitrary invertible linear transformation, but for 
simplicity, we set M z to the identity matrix. Deletion of a node is handled similarly. 

To perform an lNTERVAL p (-i/', t~, t + ) query, we first find the nodes z~ and z + storing the vertices 
Pred p (t~) and Succ p (r + ), respectively. We then successively rotate z~ with its parent until it becomes 




Figure 2: Rotation in affine composition trees. The affine function stored at a node is shown as a greek letter to its 
right (left). When rotating the pair (z\, Z2), by changing these functions as shown (right) the set of values computed 
at the leaves remains unchanged. 
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the root of the tree. Next, we do the same with z + but stop when it becomes the right child of z~ . At this 
stage, the subtree rooted at the left child of z + contains exactly all the vertices q for which q p £ [r - , t + ]. 
Thus we compose ip with the affine transformation at that node and issue the performed rotations in the 
reverse order to put the tree back in its original (balanced) position. Since z~ and z + were both within 
O(logn) steps from the root of the tree, and since performing each rotation on the tree can only increase 
the depth of a node by one, z~ is taken to the root in 0(log n) steps and this increases the depth of z + by at 
most O(logre). Thus the whole operation takes O(logn) time. 

We can augment 7(F) with additional information so that for any a G R the function E(a) = E(b\) + 
f£ F(x) dx, where b\ is the leftmost breakpoint and E(b\) is value associated with b\, can be computed in 
O(logn) time; we refer to this operation as iNTEGRATE(a). We provide the details in the Appendix. We 
summarize the above discussion: 

Theorem 3.1 A continuous piecewise-linear monotonically increasing function F with n breakpoints can 
be maintained using a data structure 7(F) such that 

1. EVALUATE and EVALUATE -1 queries can be answered in 0(log n) time, 

2. an INSERT or a DELETE operation can be performed in 0(log n) time, 

3. AFFINE and INTERVAL operations can be performed in 0(1) and 0(log n) time, respectively. 

4. INTEGRATE operation can be performed in 0(log n) time. 

One can use the above operations to compute the sum of two increasing continuous piecewise-linear 
functions F and G as follows: we first compute F(bi) for every breakpoint bi of G and insert the pair 
(bi,F(bi)) into 7(F). At this point the tree still represents 7(F) but includes all the breakpoints of G 
as degenerate breakpoints (at which the left and right derivates of F are the same). Finally, for every 
consecutive pair of breakpoints bi and bi+\ of G we apply an INTERVAL^ (ipi , bi, bj+i) operation on 7(F) 
where tpi is the affine transformation q i— > Mq+c where M = ( 1 5 ) and c = ( S) , in which Gi(x) = ax+j3 
is the linear function that interpolates between G(bi) at bi and G(bi + \) at (similar operation using 
and (i+ of G for the unbounded pieces of G must can be applied in constant time). It is easy to verify 
that after performing this series of Interval's, 7(F) turns into 7(F + G). The total running time of 
this operation is 0(|G| log \F\). Note that this runtime can be reduced to 0(|C7|(1 + log | i 7, | / log \ G\)), 
for | G | < \F\, by using an algorithm of Brown and Tarjan 00 to insert all breakpoints and then applying 
all Interval operations in a bottom up manner. Furthermore, this process can be reversed (i.e. creating 
7(F — G) without the breakpoints of G, given 7(F) and 7(G)) in the same runtime. We therefore have 
shown: 

Lemma 3.2 Given 7(F) and 7(G) for a piecewise-linear isotonic functions F and G where \G\ < \F\, 
7(F + G) or 7(F - G) can computed in 0(\G\(1 + log \F\/log\G\)) time. 



Tree sets. In our application we will be repeatedly computing the sum of two functions F and G. It will 
be too expensive to compute 7(F + G) explicitly using Lemma 3.2 therefore we represent it implicitly. 
More precisely, we use a tree set (F) = {7(F\), . . . ,7(FjS)} consisting of affine composition trees of 
monotone piecewise-linear functions F\,...,F^ to represent the function F = 2~2j=i Fj- W e perform 
several operations on F or similar to those of a single affine composition tree. 

Evaluate(x) on F takes O(klogn) time, by evaluating Ylj=i Fj( x )- Evaluate -1 (y) on F takes 
0(fclog 2 n) time using a technique of Frederickson and Johnson ifToll . 

Given the ACT 7(F ), we can convert (F) to (F + F ) in two ways: an Include(, F ) opera- 
tion sets = {7(Fi), . . . ,7(F k ),7(F )} in O(l) time. A Merge(,F ) operations sets = {7(F 1 + 
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F ),7(F 2 ), . . . , T(F fc )} in O(|F | log \Fx\) time. We can also perform the operation unInclude(, F ) 
and unMerge(, Fo) operations that reverse the respective above operations in the same runtimes. 

We can perform an Affine(, ip) where ip describes a linear transform M and a translation vector c. To 
update F by ip we update F\ by ip and for j e [2, k] update Fj by just M. This takes O(k) time. It follows 
that we can perform Interval(, Tp, r~, r + ) in O(fclogn) time, where n = \Fi\ + . . . + \Fk\. Here we 
assume that the transformation tp is such that each F{ remains monotone after the transformation. 



4 LIR on Paths 

In this section we describe an algorithm for solving the LIR problem on a path, represented as a directed 
graph P = (V, A) where V = {v%, . . . , v n } and A = {(vi, f j+i) : 1 < i < n}. A function s : V — > R is 
isotonic (on P) if s(v{) < s(vi + i), and 7-Lipschitz for some real constant 7 if s(vi + i) < s(vi) + 7 for each 
i = 1, . . . , n — 1. For the rest of this section let t : V — > R be an input function on V. For each i = 1, . . . , n, 
let Vi = {vi, . . . , Vi}, let Pj be the subpath v%,...,Vi, and let Ei and Ej, respectively, be shorthands for Ey i 
and Fp . By definition, if we let Eq = 0, then for each i > 1: 

= (x-t(^)) 2 + min S^i(y). (1) 

x— "i<y<x 



By Lemma 2.1 is convex and continuous and thus has a unique minimizer 



Lemma 4.1 For i > 1, the function E^ is given by the recurrence relation: 

{Fi l^X — 7) X ^ 5* ^ ~t~ 7 

[<_!,<-! +7] (2) 
Ei-i(x) x < s*^. 

Proof: The proof is by induction on i. E\ is clearly a single-piece quadratic function. For i > 1, since 
Ej_x is strictly convex, it is strictly decreasing on (— 00, and strictly increasing on +00). Thus 

depending on whether s*_ 1 < x — 7, s^_ 1 G [x — 7, x], or s*_ x > x, the value y G [x — 7, x] that minimizes 
Ei(y) is x — 7, s*, an d x, resp ectiv ely. ■ 



Thus by Lemmas 2.1 and 4.1 E, is strictly convex and piecewise quadratic. We call a value x that 



determines the boundary of two neighboring quadratic pieces of E\ a breakpoint of Ei. Since E\ is a simple 
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(one-piece) quadratic function, it has no breakpoints. For i > 1, the breakpoints of the function Ei consist 
of s*_ 1 and s*_ l + 7, as determined by recurrence (|2j), together with breakpoints that arise from recursive 
applications of Ei-%. Examining equation ^ reveals that all breakpoints of that are smaller than 
s*_ 1 remain breakpoints in Ei and all those larger than s*_ 1 are increased by 7 and these form all of the 
breakpoints of Ei (see Figure [3]). Thus Ei has precisely 2i — 2 breakpoints. To compute the point s* at 
which Ei(x) is minimized, it is enough to scan over these 0(i) quadratic pieces and find the unique piece 
whose minimum lies between its two ending breakpoints. 

Lemma 4.2 Given the sequence sf , . . . , s*, one can compute the j-LIR s of input function t in 0{n) time. 

Proof: For each i = 1, . . . , n, let (Tj = argmin^g^p v .\ Ei(x), i.e. U{ is the 7-LIR of Tj G where 
Tj(u) = t(v) for all j; £ 7j and in particular cr ra = s where s is the 7-LIR of t. From the definitions we get 

Ei(<Ti) = min Ei(x). 

The uniqueness of the G{ and s* implies that ai(vi) = s*. In particular, cr n (v n ) = s(v n ) = s*. Suppose that 
the numbers s(uj + x); • • ■ > s(v n ) are determined, for some 1 < i < n — 1, and we wish to determine s(uj). 
Picking s(vj) = x entails that the energy of the solution will be at least 

n 

E - *(«i)) 2 - 
j=i+l 

Thus has to be chosen to minimize on the interval [s(vi+i) — j,s(vi + i)]. If s* G — 
7, s(uj + i)], then we simply have s(uj) = s* which is the global minimum for Ei. Otherwise, one must pick 
s(vi) = s(vi + i) — 7 if s* < s(vi + \) — 7 and s(-Uj) = s(vi + \) if s* > s(vi + \). Thus each s(u,) can be found 
in 0(1) time for each i and the entire process takes 0(n) time. ■ 

One can compute the values of s*, . . . , s* in n iterations. The ith iteration computes the value s* at 
which Ei is minimized and then uses it to compute the function Ei + \ as given by (h| in 0(i) time. After 
having computed s*'s, the 7-LIR of t G can be computed in linear time. However, this gives an 0(n 2 ) 
algorithm for computing the 7-LIR of t. We now show how this running time can be reduced to 0(n log n). 

For the sake of simplicity, we assume for the rest of this paper that for each 2 < i < n, s*, i.e., the point 
minimizing Ei is none of its breakpoints, although it is not hard to relax this assumption algorithmically. 
Under this assumption, s* belongs to the interior of some interval on which Ei is quadratic. The derivative 
of this quadratic function is therefore zero at s* . In other words, if we know to which quadratic piece of Ei 
the point s* belongs, we can determine s* by setting the derivative of that piece to zero. 

Lemma 4.3 The derivative of Ei is a continuous monotonically increasing piecewise -linear function. 



Proof: Since by Lemma 4. 1 Ei is strictly convex and piecewise quadratic, its derivative is monotonically 
increasing. To prove the continuity of the derivative, it is enough to verify this at the "new" breakpoints 
of Ei, i.e. at s* and s* + 7. The continuity of the derivative can be argued inductively by observing that 
both mappings of the breakpoints from Ei-\ to Ei (identity or shifting by 7) preserve the continuity of the 
function and its derivative. 

Notice that at s*_ v both the left and right derivatives of a quadratic piece of Ei_i becomes zero. From 
the definition of Ei in ([2]), the left derivative of the function Ei — (x — tj+i) 2 at s*_ 1 agrees with the left 
derivative of at the same point and is therefore zero. On the other hand, the right derivative of the 
function Ei — (x — ii+i) 2 is zero at s*_j + 7 and agrees with its right derivative at the same point as the 
function is constant on the interval s*_ 1 + 7]. This means that the derivative of Ei — (x — U + i) 2 is 
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continuous at s*_ 1 and therefore the same holds for E{. A similar argument establishes the continuity of E% 
at s*_! + 7- ■ 
Let Fi denote the derivative of E^. Using O), we can write the following recurrence for Ff 



where 



F i+1 = 2(x-t(v i+1 )) + F i (x) (3) 

7) x>s*+7, 
Fi{x) = { xe [s|,< + 7], (4) 

X < s*. 




if we set Fq = 0. As mentioned above, s* is simply the solution of Fi(x) = 0, which, by Lemma 4.3 always 
exists and is unique. Intuitively, F is obtained from Fi by splitting it at s* , shifting the right part by 7, and 
connecting the two pieces by a horizontal edge from s* to s* + 7 (lying on the x-axis). 

In order to find s* efficiently, we use an ACT 7(F{) to represent Fi. It takes 0(log \F{\) = O(logz) 
time to compute s* = Eval uate - 1 (0) on 7(Fi). Once s* is computed, we store it in a separate array for 



back-solving through Lemma 4.2 We turn 7(Fi) into 7(Fi) by performing a sequence of Insert((s*, 0)), 
Interval x (^, a*, 00), and Insert((s*,0)) operations on where V>(<?) = q + (0); the two insert 

operations add the breakpoints at s* and s* + 7 and the interval operation shifts the portion of Fi to the 
right of s* by 7. We then turn 7(Fi) into T(-Fj + i) by performing AFFlNE(0j + i) operation on 7(Fi) where 
&+i(g) =Mg + c where M = (| °) and c = ( _ 2t (° i+l) ). 

Given ACT 'J(Fj), s* andT(Fj + i) can be computed in O(logi) time. Hence, wecan compute s*, . . . , s* 
in 0(n log n) time. By Lemma [4~2| we can conclude the following. 



Theorem 4.4 Given a path P = (V,A), a function t € M. v , and a constant 7, the j-Lipschitz isotonic 
regression of ton P can he found in 0(n log n) time. 

UPDATE operation. We define a procedure UPDATE(T(i^), t(vi±i), 7) that encapsulates the process of 
turning 7(F) into T(Fj + i) and returning s* +1 . Specifically, it performs AFFlNE(0j + i) of 7(F) to produce 
T(F i+ i), then it outputs s* +1 = Evaluate -1 ^) on 7(F i+ i), and finally a sequence of Insert((s* +1 , 0)), 
Interval^, s* +1 , 00), and Insert((s* +1 , 0)) operations on 7(F i+ i) to get 7(F i+ i). Performed on 7(F) 
where F has n breakpoints, an UPDATE takes 0(log n) time. An unUpdate(T(Fj + i), t(vi+i), 7) reverts 
the affects of an UPDATE (T(-Fj) , t(vi+i) , 7). This requires that s*^-^ is stored for the reverted version. 
Similarly, we can perform Update(, t(vi), 7) and unUpdate(, t(vi), 7) on a tree set , in 0(k log 2 n) time, 
the bottleneck coming from Evaluate" 1 (0). 

Lemma 4.5 Given 7(Fi), UPDATE(7(Fj) , t(vi), 7) and UNUPDATE(T(i^) , t(vi), 7) take 0(log n) time. 
Proof: Since Fj + i(x) = 2(x — + F(x), for a breakpoint 6^ of Fj, we get: 

b i+ i \ _ ( I ® \ ( U \ ( 



F i+ i{b l+1 ) J V 2 1 / V / V -2t(«i+i) 



where hi becomes breakpoint b 



i+l- 

We can now compute s* +l = Evaluate -1 (0) on 7(F i+ i). 

Rewriting (|4j) for a breakpoint bi + \ of Fi+i that is neither of s* +1 or s* +1 + 7, using the fact that 

i _ J + 7 h+i > 4 + i ,~ 

I bi+i < Sj+i 
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where bi + \ is the breakpoint of Pj + i that has become bi + \ in p+i, we get: 



F i+ i(b. 



■i+D 



Fi+i{bi+i) 



We next combine and rewrite equations ([5]) and (|6]) as follows: 



bi+i 
Fi+i(bi + i) 



1 



1 M F i+ i(6 i+ i) 



+ < 



%+i > s i+1 



bi+l > < + i 
bi+i < s* +1 



The new breakpoints, i.e. the pairs (s*, 0) and (s* + 7, 0), should be inserted into T(Pj). 

To perform UNUPDATE(T(Pj), t(v i), 7), we simply perform the inverse of all of these operations. 



(6) 



5 LUR on Paths 

Let P = (V, A) be an undirected path where V = {v\, . . . , v n } and A = t>i+i}, 1 < i < n} and given 
i £ M. v . For Vi £ V let Pj = (V, Aj) be a directed graph in which all edges are directed towards vi, that is, 
for j < i, (vj,Vj + i) E Ai and for j > i (vj, Vj-\) E A4. For each i = 1, . . . , n, let Tj = r(Pj, V) C IR 17 
and let Oi = argminsgPi Sy(s). If k = argminj Ey{pi), then cr K is the 7-LUR of t on P. 

We can find cr K in 0(n log 2 n) time by solving the LIR problem, then traversing the path while main- 
taining the optimal solution using Update and unUpdate. Specifically, for each i = let 
V~ = {vi, . . . , V{-i} and = {uj+i, . . . , v n }. For Pj, let Pjli, Pj+i be the functions on directed 
paths P[V^~] and P[V^ + ], respectively, as defined in Set Fi(x) = 2(x — t(vi)). Then the function 
Fi(x) = dE Vi (x) j dx can be written as Fi(x) = Fr_i(x) + F^^x) + p(x). We store Pj as the tree set 
Si = {T(F j Z 1 ),7(F^ +1 ),7(F i )}. By performing EVALUATE -1 ^) we can compute s* in 0(log 2 n) time 
(the rate limiting step), and then we can compute Ey (s*) in 0(log n) time using Integrate(s*). Assum- 
ing we have T(P i ~ 1 ) and c J(F^_ l ), we can construct T(Ff) and 7(F^_ 2 ) in O(logn) time be performing 
UPDATE(a-(i>;~i),*(vi),7) andUNUPDATE(T(P+ 1 ),t(?; i+1 ),7). Since 1 = {T(P "), T(P+), T(P)} is 
constructed in 0{n log n) time, finding k by searching all n tree sets takes 0(n log 2 n) time. 

Theorem 5.1 Given an undirected path P = (V, A) and at E MX together with a real 7 > 0, f/je ^f-LUR 
oftonP can be found in 0(n log 2 ri) time. 



6 LIR on Rooted Trees 

Let T = (V, A) be a rooted tree with root r and let for each vertex v, T v = (V V ,A V ) denote the subtree of 
T rooted at v. Similar to the case of path LIR, for each vertex v E V we use the shorthands E v = Ey v and 
E v = Ej, . Since the subtrees rooted at distinct children of a node v are disjoint, one can write an equation 
corresponding to ([T]) in the case of paths, for any vertex v of T: 

E v (x) = (x - t(v)) 2 + V min E u (y), (7) 

z — ' x—j<y<x 
u£&~ (v) 



whe re 5 (v) = {u E V \ (u, v) E A}. An argument similar to that of Lemma 4.3 together with Lemma 



2. 1 implies that for every v E V, the function E v is convex and piecewise quadratic, and its derivative F v 
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is continuous, monotonically increasing, and piecewise linear. We can prove that F v satisfies the following 
recurrence where F u is defined analogously to F in (Q: 

F v {x) = 2{x-t{v)) + P ^{x). (8) 

uG<5 — (v) 

Thus to solve the LIR problem on a tree, we post-order traverse the tree (from the leaves toward the 
root) and when processing a node v, we compute and sum up the linear functions F u for all children u 
of v and use the result in ^ to compute the function F v . We then solve F v (x) = to find s*. As in 
the case of path LIR, F v can be represented by an ACT 7(F V ). For simplicity, we assume each non-leaf 
vertex v has two children h(v) and l(v), where |^( w )| > We call h(v) heavy and l(v) light. Given 

T{Fh(v)) an d T{Fl(v))' we can compute 7(F V ) and s* with the operation UPDATE(T(F fe („) + Fh v \), t(v), 7) 
in 0(1^/(^)1(1 + log li^/^l/ log time. The merging of two functions dominates the time for the 

update. 

Theorem 6.1 Given a rooted tree T = (V, A) and a function t G MY together with a Lipschitz constant 7, 
one can find in 0{n log n) time, the ^-Lipschitz isotonic regression oft on T. 

Proof: Let T be rooted at r, it is sufficient to bound the cost of all calls to UPDATE for each vertex: 

U(T) = Ever l%)l(l + ^g |4(,)l/log l%)D- 

For each leaf vertex u consider the path to the root P u>r = (u = v\ , vi , ■ ■ ■ , v s = r ) . For each light vertex 
Vi along the path, let its value be u{v£) = 2 + \log(\F h r Vi+1 )\/\Fi( v . +1 ) =v . |)]. Also, let each heavy vertex 
Vi along the path have value v{vi) = 0. For any root vertex u, the sum of values T u = eP u r u ( v i) — 
31ogn, since \F Vi \2 v ^~ 1 < for each light vertex and \F r \ = n. 

Now we can argue that the contribution of each leaf vertex u to U(T) is at most T u . Observe that for 
!)£T, then v(l(v)) is the same for any leaf vertex u such that l(v) G P u ,r- Thus we charge u(l(v)) to v for 
each u in the subtree rooted at l(v). Since v{l{v)) > 2 + ^og{\F h ^\/\Fi^\) > 1 + log \F h ^\/\og \Fu v \\ 
(for log \Fu v \ I > log \Fu v \ |), then the charges to each v G T is greater than its contribution to U (T). Then 
for each of n leaf vertices u, we charge T u < 3 log n. Hence, U(T) < 3n log n. m 



7 LUR on Trees 

Let T = (V, A) be an unrooted tree with n vertices. Given any choice of r G V as a root, we say T( r ) is the 
tree T rooted at r. Given a function t G R v , a Lipschitz constraint 7, and a roo t r G V, the 7-LIR regression 



of i on T( r ), denoted s^ r \ can be found in O(nlogn) time using Theorem 6.1 We let £( r ) = min xg R E r (x) 

as defined in (M) and let r* = argmin rg y £( r ). The 7-Lipschitz unimodal regression (7-LUR) of t on T is 
the function > G R^. Naively, we could compute the 7-LIR in 0(n 2 log n) time by invoking Theorem 
6.1 with each vertex as the root. In this section we show this can be improved to 0(n log 3 n). 

We first choose an arbitrary vertex r G V (assume it has at most two edges in A) as an honorary root of 
T. For simplicity only, we assume that every vertex of T has degree three or less. Let the weight u(v) of a 
subtree of at v be the number of vertices below and including v in T^'K With respect to this honorary 
root, for each vertex v G V, we define p(v) to be the parent of v and h(v) (resp. l(v)) to be the heavy (resp. 
light) child of v such that u(h(v)) > uj(1(v)); f has no parent and leaf vertices (with respect to f) have no 
children. If v has only one child, label it h(v). If v = l(p(v)) we say v is a light vertex, otherwise it is a 
heavy vertex. The following two properties are consequences of oj{h{v)) > u(l(v)) for all v G T^K 
(PI) Zvev"(Kv)) = 0(nlogn). 

(P2) For u G V, the path (f, . . . , u) in contains at most [log 2 n] light vertices. 
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(a) (b) (c) 

Figure 4: Illustration of ordered with respect to honorary root f. This is redrawn, ordered with respect to structural 
root as v. The tree is redrawn with structural root as h{v) after case (a), l(y) after case (b), and p(v) after case (c). 



As a preprocessing step, we construct the 7-LIR of t for as described in Theorem 6. 1 and, for 
each light vertex v, we store a copy of 7(F V ) before merging with 7(Ff l r p r v )\). The total size of all stored 
ACTs is O(nlogn) by property (PI). We now perform an inorder traversal of T^ r \ letting each v G V 
in turn be the structural root of T so we can determine As we traverse, when v is structural root, we 
maintain a tree set (F p r v \) at p(v) with at most k = log 2 n trees, and ACTs (F h r v \) and (Fu v \) at h(v) 
and l(v), respectively. All functions F u for u = {p(v),l(v),h(v)} are defined the same as with rooted 
trees, when rooted at the structural root v, for the subtree rooted at u. Given these three data structures we 
can compute in 0(log 3 n) time by first temporarily inserting the ACTs 7(F h t v \), 7(Fu v -\), and 7(F V ) 
into (F p / V \), thus describing the function F v for a rooted tree with v as the root, and then performing an 
S W = Evaluate -1 ^) and Integrate(s (,;) ) on (F v ). Recall F v (x) = 2(x - t(v)). Thus to complete 
this algorithm, we need to show how to maintain these three data structures as we traverse . 

During the inorder traversal of T( r ) we need to consider 3 cases: letting h(v), l(v), or p(v) be the next 
structural root, as shown in Figure [4] 

(a) h(v) is the next structural root: To create (F p ^ v ^) = (F v ) we add T(F/( 1) ))) to the tree set (F p r v -\) 
through a MERGE((F p („)), 7{Fi^)) operation and perform UPDATE ((F^), t(v), 7). The ACT 
J( F l(h(v))) ls stored at l(h(v)). To create the ACT 7(F h ^ v ^) we perform unUpdate(7(F/ 1 ( w )), t(h(v)),~f), 
and then set 7(F h{h{v)) ) «- 7(F h(h{v)) - F l{h{v)) ), using 7(F l{h{v)) ). 

(b) l(v) is the next structural root: To create we perform lNCLUDE((F p ( 1 ,)), (F h ^)) and then 
Update((F p ( 1 ,)), t(v), 7). The ACTs for h(l(v)) and l(l(v)) are produced the same as above. 

(c) p(v) is the next structural root: To create 7(F V ) we set 7(F V ) <— 7(F h ^ + Fn v \) and then perform 
UPDATE('J(F 1 ,),t(v),7). To create (F p(p{v)) ) we first set (F p(j)iv)) ) <- UNUPDATE((F p(t ,)), t(p(v)), 7). 
Then if v = h(p(v)), we perform \JNMERGE((F p ^ v ^) , (Fi^ v yj)) , and 7(Fv p r v }\) is stored at 
l(p(v)). Otherwises = l(p(v)), and we perform UNlNCLUDE((# p(p(u)) ), (F h(p{v)) )) to get (F p{p{v)) ), 
and (jFh(p( v ))) is the byproduct. 

We observe that on a traversal through any vertex v G V, the number of ACTs in the tree set (F p ^) 
only increases when = u is the new structural root. Furthermore, when the traversal reverses this path 
so that p(u) is the new root where u = l(p(u)), the same ACT is removed from the tree set. Thus, when any 
vertex v G V is the structural root, we can bound the number of ACTs in the tree set (F p r v \ ) by the number 
of light vertices on the path from r to v. Thus property (P2) implies that (F p r v y) has at most log 2 n ACTs. 

We can now bound the runtime of the full traversal algorithm described above. Each vertex v G V 
is visited as the structural root 0(1) time and each visit results in O(l) Update, unUpdate, Include, 
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unInclude, Merge, and unMerge operations. Performing a single Update or unUpdate on a tree 
set takes at most 0(fclog 2 n) = 0(log 3 n) time, and each INCLUDE and UNlNCLUDE operation takes 0(1) 
time, so the total time is 0(n log 3 n). We now need to bound the global costs of all MERGE and unMerge 
operations. In proving Theorem |6.1| we showed that the cost of performing a MERGE of all light vertices 
into heavy vertices as we build (F r ~) takes 0(ralog 2 n). Since at each vertex v € V we only MERGE and 
unMerge light vertices, the total time of all Merge and unMerge operations is again 0(n log 2 n). 
Thus we can construct £W for each possible structural root, and using r* = arg min„ g y £W as the 



optimal root we can invoke Theorem 6. 1 to construct the 7-LIR of t on T^ r *\ 



Theorem 7.1 Given an unrooted tree T = (V, A) and a function t E R v together with a Lipschitz constraint 
7, we can find in 0(n log 3 n) time the ^-Lipschitz unimodal regression oft on T. 



8 Conclusion 

In this paper we defined the Lipschitz isotonic/unimodal regression problems on DAGs and provided effi- 
cient algorithms for solving it on the special cases where the graph is a path or a tree. Our ACT data structure 
has proven to be a powerful tool in our approach. 

Our algorithms can be generalized in a number of ways, as listed below, without affecting their asymp- 
totic running times. 

• One can specify different Lipschitz value 7 for each Lipschitz constraint, thus writing the Lipschitz 
constraint involving Sj and Si + \ as 8i+i < S{ + 7$. 

• Isotonicity constraints can be regarded as (backward) Lipschitz constraints for which the Lipschitz 
constant is zero. One can allow this constant to be chosen arbitrarily and indeed permit different 
constants for different isotonicity constraints. In other words, the constraint Sj < Sj + i can be replaced 
with s i+ i > Si + 5i. 

• The £2 norm || • || can be replaced with a "weighted" version || • || \ for any A : V — ► M by defining for 
a function s : V — > M: 

\\s\\l = J2 X (v)s(v) 2 . 

The most important open problem is solving LIR problem on general DAGs. The known algorithm for 
solving IR on DAGs runs in 0(n 4 ) and does not lend itself to the Lipschitz generalization. The difficulty 
here is to compute the function Ei (or in fact Fj) when the zth vertex is processed and comes from the 
fact that the function, j = 1, . . . , k, where i\, . . . , are vertices with incoming edges to i, are not 
independent. Independently, it is an important question whether IR can be solved more efficiently. In 
particular, the special case of the problem where the given DAG is planar has important applications in 
terrain simplification. 
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A Integrate Operation 



Lemma A. 1 We can modify an ACT datastructure for 7(F) with \F\ = n and F(x) = dE(x) / dx, so for 
any x G Rwe can calculate E(x) = iNTEGRATE(x) in 0(log n) time, and so the cost of maintaining 7(F) 
is asymptotically equivalent to without this feature. 

Proof: For any a G R, we calculate iNTEGRATE(a) = E(a) by adding E(b\) to J* fe a F(x) dx, where b\ 
is the leftmost breakpoint in F. Maintaining E(b\) is easy and is explained first. Calculating f b F(x) dx 
requires storing partial integrals and other information at each node of 7(F) and is more involved. 

We explain how to maintain the value E(b\) in the context of our 7-LIR problems in order to provide 
more intuition, but a similar technique can be used for more general ACTs. We maintain values ci, C2, C3 
such that ax 2 + C2X + 03 = Yli=i( x ~ ^*) 2 f° r tne & = 0(n) values of ij. This is done by adding 1, — 2ti, 
and tf to c\, C2, and C3 every time a new t{ is processed. Or if two trees are merged, this can be updated in 
O(l) time. The significance of this equation is that E(b\) = c\b\ + C261 + 03 because the isotonic condition 
is enforced for all edges in A, forcing all values Sj in the 7-LIR to be equal. 

When we construct 7(F) we keep the following extra information at each node: each node z whose 
subtree spans breakpoints bi to bj maintains the integral I(z) = J b J ^ z 1 (F(x) — F(bi)) dx (this is done 
for technical reasons; subtracting F(bi) ensures F(x) — F(bi) > for x G [bi, bj] since F is isotonic), the 
width W(z) = ^(bj)-®- 1 ^), and the height H(z) = ^- 1 (F(b j ))-^- 1 (F(b i )). Recall that$ 2 is the 
composition of affine transformations stored at the nodes along the path from the root to z. We build these 
values (I(z), W(z), and H(z) for each z G 7(F)) from the bottom of 7(F) up, so we apply all appropriate 
affine transformations cj) z . for decedents zj of z, but not yet those of ancestors of z. In particular, if node z 
has left child z\ and right child z k , then we can calculate 

/(*) = cp zi (I( Zl )) + <t> tk (I(z k )) + <p Zl (H( Zl )) ■ 4>z k (W(z k )) (see Figure 
H(z) = cf )zi (H(z l )) + ( t }zk (H(z k )), and 
W(z) = ct> Zt (W(z l )) + cl )zk (W(z k )). 

When a new breakpoint b is inserted or removed there are O(logn) nodes that are either rotated or are on 
the path from the node storing b to the root. For each affected node z we recalculate I(z), W(z) and H(z) 
in a bottom up manner using z's children once they have been updated themselves. Likewise, tree rotations 
require that 0(1) nodes are updated. So to demonstrate that the maintenance cost does not increase, we just 
need to show how to update I(z), H(z) and W(z) under an affine transformation 4> z (z) in 0(1) time. 

For affine transformation <f) z applied at node z (spanning breakpoints bi through bj), we can decompose 
it into three components: translation, scaling, and rotation. The translation does not affect I(z), W(z), 
or H(z). The scaling can be decomposed into scaling of width w and of height h. We can then update 
W(z) <— w • W(z), H(z) = h • H(z), and I(z) = w ■ h ■ I(z). To handle rotation, since translation does 
not affect I(z), W(z), or H(z), we consider the data translated by ^~ 1 ((-b i , 0)) so that F(bi) - 

F(bi))) = (0,0) is the origin and ^((bj, F(bj) - F(bi))) = (x,y) = (W(z),H(z)). Then we consider 
rotations about the origin, as illustrated in Figure [5] After a rotation by a positive angle 6, the point (x, y) 
becomes (x',y ! ) = (xcosO — ys'm8,xsiTi9 + ycos9). This updates W(z) <— x' and H(z) <— y' . 
We can update I(z) by subtracting the area A_ that is no longer under the integral y • (ytan#)/2 and 
adding the area A + newly under the integral x' ■ (x'tan^)/2 as illustrated in Figure [5] Thus I(z) <— 
I(v) + x 1 ■ (x'tan#)/2 — y ■ (ytan#)/2. A similar operation exists for a negative angle 9, using the fact 
that F must remain monotone. As desired all update steps take 0(1) time. 

Now for a value a G R we can calculate iNTEGRATE(a) = E(a) as follows. Let z a G 7(F) be the 
node and b be the breakpoint associated with PRED^(a). We consider the path Z = (z a , . . . , zi, Zq) from z a 
to the root z$ of 7(F), and will inductively build an integral along Z. We calculate I Zo = J fo a 4> ZQ (F(x) — 
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Figure 5: Left: Illustration of the change in integral caused by a rotate of 9 degrees counterclockwise. Right: Calcula- 
tion of I z from children Zk and zi of z. 

F(bi)) dx similar to the way we calculated I(zq), but just over the integral [61, a]. As a base case, let 
I Za = &~^(F(x) — F(b))dx (which is easy to calculate because F is linear in this range) and let 
W Za = $~ 1 (a) — f»J (6) describe the integral and width associated with the range [b, a] and z a . For a node 
z S Z that spans breakpoints bi through bj, we calculate an integral I z = f^ 1 ^~ 1 (F(x) — F(b{)) dx and 
width W z = < I ) J 1 (a) — $J 1 (6j) using its two children: z\ and z^. Let be the right child which lies in Z 
and for which we had inductively calculated I Zk andW Zk . Thenlet/^ = <fi zl (I(zi))+<f) Zk (I Zk )+<fi zl (H(zi))- 
(f>z k (W Zk ) (as illustrated in Figure[5]> and let W z = <j) Zl {W{zi)) + <f> Zk (W Zk ). Once (j> Zo (I Zo ) and (f) Zo (W Zo ) 
are calculated at the root, we need to adjust for the F(bi) term in the integral. We add 4> ZQ (W Z0 ) ■ F(bi) 
to (f> Zo (I(z )) = f£(F(x) - F(bi)) dx to get J & a F{x) dx. So finally we set iNTEGRATE(a) = E{b\) + 
4>z (Iz ) + 4>z (W Zo ) ■ F{b\). Since the path Z is of length at most O(logn), E(a) can be calculated in 
0(log n) time. ■ 
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